An integrative systems biology approach to overcome venetoclax resistance in acute myeloid leukemia

The over-expression of the Bcl-2 protein is a common feature of many solid cancers and hematological malignancies, and it is typically associated with poor prognosis and resistance to chemotherapy. Bcl-2-specific inhibitors, such as venetoclax, have recently been approved for the treatment of chronic lymphocytic leukemia and small lymphocytic lymphoma, and they are showing promise in clinical trials as a targeted therapy for patients with relapsed or refractory acute myeloid leukemia (AML). However, successful treatment of AML with Bcl-2-specific inhibitors is often followed by the rapid development of drug resistance. An emerging paradigm for overcoming drug resistance in cancer treatment is through the targeting of mitochondrial energetics and metabolism. In AML in particular, it was recently observed that inhibition of mitochondrial translation via administration of the antibiotic tedizolid significantly affects mitochondrial bioenergetics, activating the integrated stress response (ISR) and subsequently sensitizing drug-resistant AML cells to venetoclax. Here we develop an integrative systems biology approach to acquire a deeper understanding of the molecular mechanisms behind this process, and in particular, of the specific role of the ISR in the commitment of cells to apoptosis. Our multi-scale mathematical model couples the ISR to the intrinsic apoptosis pathway in venetoclax-resistant AML cells, includes the metabolic effects of treatment, and integrates RNA, protein level, and cellular viability data. Using the mathematical model, we identify the dominant mechanisms by which ISR activation helps to overcome venetoclax resistance, and we study the temporal sequencing of combination treatment to determine the most efficient and robust combination treatment protocol.


A.1 Drug uptake and regulation of transcription factors
The first two modules, cellular drug uptake and the drug effects on transcription factor expression, were modeled as described in Ref. [1]. Here we give a brief overview of the model and refer the reader to Ref. [1] for more details. In this approach, the in vitro cellular uptake and decay of drug D = {A, T } for venetoclax or tedizolid, respectively, was modelled as a system of two coupled ordinary differential equations (ODEs) [1]: The mean-field model with temporal delay [1] was utilized in this work to model the expression of c-Myc and Chop in the presence of drug administration. Following this convention, we will refer to the transcription factors c-Myc and Chop as "target proteins", and their concentration as [T j ]. In this approach, the temporal evolution of [T j ] are modeled by treating the drugs as transcriptional regulators, assuming the quasi-steady state for mRNA expression levels [2]. The monotherapy with either  or Tedizolid perturbs the regulatory networks of the target proteins, resulting in four distinct possible effects on the target protein expression, which are: (1) increase in production rate, (2) decrease in production rate, (3) increase in stability, and (4) decrease in stability. Rather than distinguishing between the underlying mechanisms leading to each drug effect, an appropriate (mean) contribution is included in the model to capture the overall regulatory behaviour. For a single therapeutic agent, the time evolution of [T j ] is given by the following first-order ODE [1]: In the above equation j = 1, 2 corresponds to c-Myc or Chop, respectively, subscript m labels the drug (m = 1, 2 corresponds to venetoclax or tedizolid, respectively), k Tj ,0 is the background production rate for target protein T j in the untreated case, k Tj ,1 is a multiplier that dictates the maximum production rate for T j in the presence of drug, and δ Tj ,0 , δ Tj ,min , and δ Tj ,max are, respectively, the nominal, minimum and maximum decay rates for T j . We note that k Tj ,0 , k Tj ,1 , δ Tj ,0 , δ Tj ,min , and δ Tj ,max are biologically-constrained parameters that do not depend on the drug concentration. The drug concentration is contained in the terms H This is taken to be the initial condition for the treated system in numerical simulations. Following Ref. [1], the effects of combination therapy on the time evolution of [T j ] is described by the ODE: j,ℓ and H (2) j,ℓ , which are defined by: and H (2) j,ℓ = 1 − H (1) j,ℓ = 1 1 + C j,ℓ , (A. 8) where C j,ℓ = 2 m=1 w j,m,ℓ [X j,m,ℓ ] n j,m,ℓ + q j,ℓ 2 m=1 w j,m,ℓ [X j,m,ℓ ] n j,m,ℓ . (A.9) It is clear that H (1) j,ℓ and H (2) j,ℓ are Hill functions of the argument C j,ℓ , which as Eq. (A.9) illustrates, contains all instances of the drug concentrations, [X j,m,ℓ ]. Following the work of Ref. [5], the weights w j,m,ℓ , where w j,m,ℓ ≥ 0 ∀ {j, m, ℓ}, were introduced to incorporate antagonistic drug effects. The weights distinguish the relative importance of the individual drug effect to the combined effect during combination therapy: taking w j,m,ℓ = 1 corresponds to each drug contributing equally to the combined effect (the default value), while taking w j,m,ℓ ≈ 0 captures cases where combined antagonistic drug effects nullify or dampen a pathway that plays a key role in the individual response to one or both drugs. Importantly, taking w j,m,ℓ = 1 in Eq. (A.9) gives the familiar expression for protein regulation by two transcription factors, where cooperativity between transcription factor binding is described by the constant factor q j,ℓ [2,5]. This cooperativity between two or more drugs might arise when the protein networks affected by each drug overlap or have synergistic physiological functions [1].
Since the mean-field approach does not model individual molecular pathways, temporal delays can be incorporated into the drug effects to account for regulatory interactions that are downstream of several other protein interactions. For example, if the administration of a drug affects the expression of a protein that is upstream of a regulatory pathway, there may be a series of interactions, including transcription and translation events and post-translational modifications, that must occur before the regulatory protein is affected. In this scenario, there should be a time delay between drug administration and the corresponding effect on the target protein. To incorporate temporal delays into the drug effects, we make the replacement [X j,m,ℓ ](t) = [D m,i ] t − τ j,m,ℓ /K j,m,ℓ , where τ j,m,ℓ is the time delay. The kinetic parameter values that were previously determined for the MOLM-13 R2 cell line using the mean-field approach with temporal delay are presented in Table D.

A.2 Regulation of Bcl-2 family proteins by transcription factors
In Table A, we present the direct and indirect individual regulatory effects of the transcription factors c-Myc and Chop on the expression of the Bcl-2 family protein levels that were included in the model. Square brackets denote the concentration of the corresponding protein, t denotes the time, and d[X] dt denotes the time-derivative of the concentration of protein X. In addition, k 0,X , k X , n i , and K i are constants. The constant k 0,X represents the background (basal) production rate for protein X, and the maximal production rate is given by k 0,X + k X . We use Hill functions [3,4] to model the regulatory effects. Furthermore, we make the quasi-steady state approximation for mRNA levels so that protein production can be treated as a reduced one-step process [2]. The shape and steepness of the Hill function for regulatory effect i is determined by the value of the exponent n i , which is referred to as the Hill coefficient due to its origin in the context of ligand binding to macro-molecules [3]. The constant K i is the half-saturating concentration for the transcription factor in the regulatory event i. Further details about the regulatory interactions are described in Section 2.1.3 of the manuscript for c-Myc and Section 2.1.4 for Chop.
As indicated in Table A, each of the Bcl-2 protein family members is regulated by both c-Myc and Chop. Furthermore, due to the bistable response of c-Myc [8,21], some of the Bcl-2 protein family members experience two different regulatory mechanisms due to c-Myc alone, depending on the c-Myc concentration. When multiple transcription factors regulate the expression of a gene, their effects must be considered simultaneously. To model the combined regulatory effects of c-Myc and Chop, including the bistable response of c-Myc [8,21], we follow the work of Refs. [2,5]. The only exception to this is for the combined regulatory effects of Bim since the inhibitory effect of c-Myc on Bim expression is the prevention of Bim accumulation in the mitochondria, see Section 2.1.3 of the manuscriptand Refs. [6,9,10]. Thus, even if Chop concentration is sufficient to induce Bim expression, it will not accumulate in the mitochondria to play a role in the apoptosis pathway, and this is taken into consideration in the ODE for Bim, c.f. Eqs. (A.20) and (A.21) and their following descriptions, which are presented in Section A.7, after first discussing the additional components of the model. We note that in implementing the convention in Refs. [2,5] for the combined regulatory effects of multiple transcription factors, we make the following assumptions. (i) Each Bcl-2 family member has a promoter with at least two non-overlapping operator regions, and each transcription factor can bind to its own operator region to regulate the expression of a target gene. (ii) The rate of gene expression contains a constant background (basal) production term that is independent of the concentration of the transcription factors, denoted by k 0,X for protein X. (iii) Following the hypothesized bistable response of c-Myc [8,21], c-Myc is an activator of gene expression for Mcl-1 and Bcl-2 at lower concentrations and a repressor at higher concentrations; similarly, c-Myc is a repressor of gene expression for Bim at lower concentrations and an activator of Bim, Bax, and Bak expression at higher concentrations. This amounts to assuming that K 2 < K 4 and K 3 < K 5 . We do not impose additional constraints on the half-saturation constants or define a hard cut-off that segregates the hypothesized cancer zone and apoptosis zone. (iv) If a repressor is bound to its operator region, activated gene expression is fully blocked, except for the basal expression level.
Finally, we point out that our model neglects the inhibitory effects of ISR activation on protein translation rates. We justify this assumption by noting that, under experimental conditions, EIF2α is dephosphorylated by GADD34 within hours of ISR activation [22], thus global translational  [19] inhibition occurs over a short time scale. In comparison, the experiments considered in this work span several days, and we therefore expect the inhibitory effects of EIF2-α phosphorylation on the expression of the proteins in our model to be minimal over the long-term. As discussed in Section 2.1.4 of the manuscript, it has recently been proposed that there are two nodes for translational inhibition in the ISR. In particular, in a delayed response to prolonged conditions of cellular stress, transcriptional and post-translational regulatory effects shift the cellular response toward cap-independent translation [22]. Importantly, this results in the selective translation of stress and apoptosis proteins, in addition to proteins required for regulatory cell processes, such as c-Myc [22][23][24]. While Mcl-1 was seen to be downregulated when cap-dependent translation was inhibited in human leukemia cells [25], we note that this effect is implicitly captured in the model by the inhibitory effect of Chop on Mcl-1 expression. Importantly, these simplifications and assumptions are consistent with the long-term effects of ISR activation and also prevent the unnecessary introduction of additional kinetic parameters into the model.

A.3 Bcl-2 family protein interactions
We limit the scope of the model to the Bcl-2 family proteins for which treated experimental expression levels were available. Furthermore, as noted in Section 2.3 of the manuscript, we excluded certain Bcl-2 family proteins from the model based on redundancy in protein function and also based on their lower expression in the resistant cell line, indicating a limited role in the development of resistance to venetoclax. This enabled a significant decrease in the number of kinetic parameters in the model, while still retaining sufficient detail to capture the main interactions in the system. Specifically, the model includes the expression of, and interactions between, Bcl-2, Mcl-1, Bim, Bax, and Bak. The interactions, including the reaction equations and corresponding rate laws, are presented in Table B. We note that protein complexes are indicated by the notation [X:Y ], and k i , k −i , and κ i are constants which denote, respectively, association rates, dissociation rates, and activation rates. Table B: Bcl-2 family protein-protein interactions comprising the apoptosis pathway in the model.

Description
Reaction equation

Rate laws Reference
Mcl-1 binds to and dissociates from Bim [ [26][27][28][29][30] Bcl-2 binds to and dissociates from Bim [ [26][27][28][29][30] Bim binds to and activates Bax Bim binds to and activates Bak Mcl-1 binds to and inhibits activated Bax [26,27,29,30,34] Mcl-1 binds to and inhibits activated Bak [26,27,29,30,35,36] Bcl-2 binds to and inhibits activated Bax [26,27,29,30,37] Bcl-2 binds to and inhibits activated Bak [26,27,29,30,35] A.4 Caspase activation and feedback As discussed in Section 2.1.1 of the manuscript, the activated pro-apoptosis effector proteins Bax and Bak form homo-and heterodimers and higher order oligomers which embed into the mitochondrial outer membrane (MOM). This promotes mitochondrial outer membrane permeabilization (MOMP) and the release of cytochrome C into the cytosol, which is followed by the formation of the apoptosome, the activation of several Caspases and subsequent negative feedback on the prosurvival Bcl-2 family proteins, and ultimately cell death. To model these processes, we use a simplistic approach that enables a reduction in the number of equations and kinetic parameters in the mathematical model. Specifically, we do not directly model the formation of active Bax and Bak homo-and heterodimers and oligomers, the process of MOMP, or the Caspase cascade. Rather, as indicated in Table C, we take the level of active Caspase-3 to be indicative of the onset of apoptosis since Caspase-3 is known to be an executioner of apoptosis. Furthermore, we make the simplifying assumption that the level of active Caspase-3 depends directly on the concentrations of active Bax and Bak in the system. Finally, we include the negative feedback of Caspase-3 on the pro-survival Bcl-2 family proteins, which amplifies the apoptotic signalling at the mitochondria.
Activated Bak promotes MOMP and caspase activation Activated caspase 3 cleaves and inhibits Bcl-2 A.5 Interaction of venetoclax with Bcl-2 family proteins Venetoclax (ABT-199) is a potent and specific inhibitor of the Bcl-2 protein. It works by sequestering and binding to Bcl-2 and forming a protein-drug complex with a high binding affinity. This prevents Bcl-2 from carrying out its anti-apoptotic functions, specifically, from binding to the BH3-only proteins and the pro-apoptosis effector proteins. We model this interaction as follows: where k 13 and k −13 are, respectively, the association and dissociation rates. The corresponding rate laws for this interaction are given by: -199] denotes the intracellular concentration of venetoclax that is not bound to Bcl-2.
In addition to the regulatory effects induced by c-Myc and Chop that resulted from the administration of venetoclax, the protein expression data indicated that venetoclax treatment causes a significant decrease in the expression of Bax. Importantly, this downregulation in Bax expression could not be explained by the direct and indirect regulatory effects of c-Myc and Chop outlined in Table A. This is because on day 3 of ABT-199 treatment, c-Myc expression has increased, see Figure 3 of the manuscript. Since c-Myc upregulates Bax expression, we would expect to see a resulting increase in Bax levels during treatment. On the other hand, if a decrease in Chop expression were the cause of the drop in Bax levels on day 3 of ABT-199 treatment, then the same mechanism would cause a decrease in Bax levels with tedizolid treatment as well, which is inconsistent with experimental measurements, see Figure 3 of the manuscript. Consequently, there is another mechanism that is leading to the marked decrease in Bax expression during ABT-199 treatment.
As discussed in Section 2.1.2, activation of the PI3K/AKT pathway has been commonly observed in cells after both short and long-term exposure to venetoclax treatment [49]. Furthermore, it has been previously reported [50] that AKT phosphorylates Bax at residue S184, which impedes its ability oligomerize and function as a pro-apoptosis protein and also significantly decreases its protein stability [51,52]. Since these changes provided biological rationale for the downregulation of Bax expression with venetoclax treatment, we included an inhibitory term in the model which depends on the concentration of venetoclax. We point out that this inhibitory term was taken to have the form of a Hill-function and was applied to the production rate of Bax to capture the collective changes of Bax phosphorylation: where K 13 and n 13 are constants, and [ABT-199 total ] denotes the total intracellular concentration of venetoclax. Importantly, we include this inhibitory interaction when combining the regulatory effects of the transcription factors on Bax expression in Section A.7. We note that, alternatively, the effects of Bax phosphorylation by AKT could potentially be captured by including a term in the model which decreases the half-life of the Bax protein in a venetoclax-concentration dependent manner. The limitation of this approach is that an additional kinetic parameter must be introduced into the model, namely, the half-life of the phosphorylated protein. Furthermore, this term alone may not collectively capture all of the changes that occur as a result of Bax phosphorylation at S184.

A.6 Cellular proliferation and viability
We use a simplistic approach to model cellular proliferation and cellular death that takes into account the main mechanisms discussed in Sections 2.1.1 -2.1.4 of the manuscript. Specifically, we assume that the time evolution of the number of live cells in the model can be described by the equation: is the rate of cellular proliferation and β ≡ β [Casp-3] is the rate of cell death. The rate of cellular proliferation depends on several variables to account for the major metabolic effects of the proteins depicted in the pathway in Figure 2 of the manuscript. The functional form for α is taken to be: (A.14) where λ 0 , λ Bcl-2 , λ c-Myc , K 14 , K 15 , K 16 , n 14 , n 15 , and n 16 are constants. In Eq. (A.14), the first term describes a constant background production that is independent of the protein expression levels in the model, and the second term describes the stimulatory effect of Bcl-2 on cellular (mitochondrial) respiration [49,[53][54][55][56][57][58][59], see Section 2.1.1. This term implicitly takes into account the inhibitory effects of venetoclax treatment on cellular respiration since administration of venetoclax causes a significant drop in free Bcl-2 levels. In addition, the cellular proliferation depends on c-Myc concentration to account for the stimulatory effects of c-Myc on glycolysis and cellular respiration [60][61][62][63][64][65][66][67][68][69], see Section 2.1.3. Finally, we include an inhibitory term that depends on the intracellular tedizolid concentration, [Ted], to account for the repression of cellular respiration previously observed, see Refs. [59,70] and Section 2.1.4. The cellular death rate is taken to depend upon the level of active Caspase-3 [29,30,38,39,[71][72][73][74][75][76][77][78][79][80], with the following form: where Ω Casp-3 , K 17 and n 17 are constants, and [Casp-3* frac ] is the fraction of total Caspase-3 in the system that is active, . While technically not a concentration, we have maintained the square brackets on the quantity for notational convenience. Using this model, the time evolution of the number of dead cells in the system is described by: 16) and the cellular viability is calculated as follows: For consistency with experiments, see Section 2.2, in numerical simulations we initialized the number of live cells to 50,000. We initialized the number of dead cells in numerical simulations such that the initial cellular viability was equal to the average of the experimentally measured untreated cell viability. This amounted to initializing the number of dead cells in numerical simulations to 813, corresponding to an untreated cellular viability of 98.4%. Proceeding this way maintained an approximately constant cellular viability in the untreated simulations, while still capturing the net cellular proliferation observed in experiments. In the case of treatment simulations, drug-induced fluctuations in the protein levels led to time-dependent changes in the cellular proliferation and death rates, and consequently, also in the simulated cellular viabilities.

A.7 System of coupled ODEs for the apoptosis proteins
Now that we have mapped out the individual interactions in the pathway presented in Figure 2, we combine all of the individual rate laws for each apoptosis protein. This leads to a system of coupled ODEs that describe the time evolution of the concentration of each protein species in the network. Below we describe the set of equations.
1. The equation describing the rate of change of free Mcl-1 concentration is given by: The first term is the basal production rate and the second term describes the stimulatory regulatory effect of c-Myc and the inhibitory effect of Chop. The next six terms describe binding and dissociation events with other Bcl-2 family proteins, the seventh term describes cleavage by active Caspase-3, and the last term describes natural protein decay.
The first term is the basal production rate and the second term describes the stimulatory regulatory effect of c-Myc at lower concentrations, the inhibitory regulatory effect of c-Myc at higher concentrations (K 2 < K 4 ), and the inhibitory effect of Chop. The next six terms describe binding and dissociation events with other Bcl-2 family proteins, the seventh and eighth terms describe binding and dissociation with venetoclax, the ninth term describes cleavage by active Caspase-3, and the last term describes natural protein decay. The first term is the basal production rate and the second term describes the inhibitory regulatory effect of c-Myc at lower concentrations, the stimulatory regulatory effect of c-Myc at higher concentrations (K 3 < K 5 ), and the stimulatory effect of Chop. We note that the production terms give dominance to the regulatory effects of c-Myc on the expression of Bim. In particular, the terms in square parentheses capture protein production by c-Myc when [c-Myc] > K 5 , in which case the stimulatory Hill function will saturate, the maximal production rate k 0,Bim + k Bim is achieved, and production by Chop is redundant. For lower concentrations of c-Myc, K 2 < [c-Myc] < K 5 , the Hill function describing production by c-Myc will be close to zero, but so will the Hill function describing the inhibitory effect of c-Myc, which will ultimately lead to a low rate of production for Bim. In this case, the production of Bim by Chop is also blocked. We chose this form for the production rate of Bim due to the indirect nature of the inhibitory effect caused by c-Myc. In particular, in this limit, c-Myc prevents Bim translocation to the mitochondria. Thus, even if Chop concentration is sufficient to induce Bim expression, it will not accumulate in the mitochondria to play a role in the apoptosis pathway. Finally, at lower concentrations of c-Myc, [c-Myc] < K 2 , the Hill function describing production by c-Myc will be close to zero, but the Hill function describing inhibition by c-Myc will be close to one, i.e. the inhibitory effect will be negligible. In this case, expression by Chop will dominate the production of Bim and if Chop expression is high enough, the maximal production rate will be achieved.
The next four terms describe binding and dissociation events with the pro-survival Bcl-2 family proteins, and the following six terms describe binding, dissociation, and activation of the pro-apoptosis effector proteins. The last term describes natural protein decay.
4. The equation describing the rate of change of free Bax concentration is: The first term is the basal production rate and the second term describes the stimulatory regulatory effect of c-Myc at high concentrations, the stimulatory regulatory effect of Chop, and the inhibitory effect of venetoclax. The next two terms describe binding and dissociation events with Bim, and the last term describes natural protein decay. The first term is the basal production rate and the second term describes the stimulatory regulatory effect of c-Myc at high concentrations, and the stimulatory regulatory effect of Chop. The next two terms describe binding and dissociation events with Bim, and the last term describes natural protein decay. The first term describes the rate of activation by pro-apoptosis activator protein Bim, the next four terms describe binding and dissociation events with pro-survival Bcl-2 family proteins, and the last term describes natural protein decay. The first term describes the rate of activation by pro-apoptosis activator protein Bim, the next four terms describe binding and dissociation events with pro-survival Bcl-2 family proteins, and the last term describes natural protein decay.

The equation describing the time evolution of (inactive) Caspase-3 concentration is:
The first term is the production rate of Caspase-3, the next two terms describe the activation of Caspase-3 by the activated pro-apoptosis effector proteins, and the last term describes natural protein decay. The first term describes the cleavage of Mcl-1 by active Caspase-3 and the second term describes natural protein decay.
11. The equation describing the time evolution of inactive (cleaved) Bcl-2, denoted Bcl-2 − , concentration is: The first term describes the cleavage of Bcl-2 by active Caspase-3 and the second term describes natural protein decay.
12. The equations describing the time evolution of the concentration of protein complexes in the system are given by: In the above equations, the first term describes complex formation, the second term describes complex dissociation, and the last term describes natural protein decay. In Eqs. (A.37) In the above equation, the first term describes the formation of the protein-drug complex, the second term describes the complex dissociation, and the last term describes natural decay, which is assumed to occur at the same rate as the free intracellular venetoclax concentration, We point out that each combined ODE contains a constant decay rate. This decay is assumed to be the result of natural spontaneous degradation as well as dilution in concentration caused by cell growth [2]. We also note that since each protein in the apoptosis pathway exists in several forms, to obtain the total concentration of each protein, we must sum up the concentration for all of its species. For example, the time-dependent total concentration of Mcl-1 is given by: Finally, we note that the full system of coupled ODEs that describes the time evolution of the coupled intrinsic apoptosis-ISR pathway, Figure 2, is given by Eqs. (A.9)-(A.28), along with the equations for cellular drug uptake and the regulation of transcription factors described in Section A.1, and the equations for cellular proliferation and viability described in Section A.6.

B Kinetic parameter values
The values of the kinetic parameters for the cellular drug uptake model and the regulation of the transcription factors c-Myc and Chop, see Section A.1, were previously determined for the MOLM-13 R2 AML cell line [1] and are presented in Table D. Several of the remaining kinetic parameter values were determined directly from the literature. In Table E, we present the values of the free protein decay rates and protein complex decay rates that were obtained this way. hour −1 [83,89] In Table F, we present the forward and reverse rate constants for the formation of protein complexes in the apoptosis pathway. In addition, we present the corresponding dissociation constants obtained from the literature, as well as the activation rate constants for the pro-apoptosis effector proteins. Table F: Binding, dissociation, and activation rates for the apoptosis proteins, measured for various cell lines and experimental conditions, determined from the literature. Seconds is abbreviated by "s" and nanomolar concentration is abbreviated by "nM". Following Ref. [83], the dissociation constants K D and backward constants were obtained from the literature where available and the forward constant was calculated as k forward = k backward /K D .  [83,89,91] The remaining parameters were determined by fitting to the experimental data, where the fitting procedure, as well as the procedure for initializing and performing numerical simulations, is described in Section 2.4 of the manuscript. The nominal parameter set determined this way, along with the corresponding parameter search ranges, are presented in Table G. We approximated the search ranges for the initial levels of Bcl-2 family proteins based on ranges that were observed via quantitative Western blotting for several cancer cell lines and patient tumor samples in Ref. [83]. We used a slightly larger range for Mcl-1 levels than were previously observed since our protein data indicated that the upregulation of Mcl-1 was a significant contributor to venetoclax resistance in the MOLM-13 R2 cell line (it was approximately 40 times higher in the resistant versus parental cell line). We used a large range for the half-saturation constants to account for the possibility that some interactions that have been reported in the literature may not be the dominant regulatory interactions in the system. We note that the kinetic parameter values for the equations describing cellular drug uptake and the regulation of transcription factors, see Section A.1, were reported in Supplemental Table 2 of Ref. [1].

Parameter
As explained in Section 2.4 of the manuscript, the equilibrium simulations were initialized by first turning off all protein-protein interactions, except for the production and decay terms [92]. Thus the initial values of all protein complexes, active pro-apoptosis effector proteins and Caspase levels, and cleaved pro-survival protein levels were taken to be zero. Consequently, the free protein levels, whose time evolution is described by Eqs. (A.18)-(A.22) and (A.25), were initialized to the total initial protein levels reported in Table G. The transcription factor-dependent production rates were calculated from the background production rate, Hill functions, and decay rates, see Table G, for each protein in the apoptosis pathway by making the quasi-steady state [2] assumption. For Mcl-1, for example, using Eq. (A.18), this leads to: for the values presented in Table G. For the other Bcl-2 family proteins, we find k Bcl-2 = 2.52 nM/min, k Bim = 0.02195 nM/min, k Bax = 0.5002 nM/min, and k Bak = 0.8604 nM/min. For Caspase-3, the production rate was set to balance the decay rate since there was no assumed transcription factordependent production of this protein.
Lastly, in Table G, we have reported individual search ranges for each rate that contributes to the total cellular proliferation. However, we note that when searching for these parameter values, it was the total doubling rate in the untreated system (which includes multiplication by Hill functions, see Eqs. (A.14)-(A.15)) that was constrained to fall within these previously reported biologically-relevant ranges [93]. In particular, given the value of the initial protein levels, halfsaturation constants, Hill coefficients, and proliferation and death rates in Table G, the contribution to cellular proliferation from Bcl-2 is at most 1.09 × 10 −3 hour −1 and from c-Myc it is 3.39 × 10 −5 hour −1 in the untreated system. This leads to a maximum rate of cellular proliferation in the system of 4.5 × 10 −2 ≈ ln(2)/15.3 hour −1 , i.e. a cell division time of about 15.3 hours, which corresponds to twice the minimum doubling time observed previously in acute leukemia cells [93]. Moreover, to achieve the level of active Caspase-3 that was observed experimentally in the treated system, we imposed a maximum of 0.15 for the active Caspase-3 fraction in the untreated system, leading to an upper limit on the death rate in untreated cells of 7.62 × 10 −3 ≈ ln(2)/91.0 hour −1 . Taking into account cellular proliferation and cell death, this corresponds to a doubling time of approximately 18.4 hour −1 in the untreated system, which agrees with the experimental data used in this study and falls well within previously reported biologically-relevant ranges.

C Sensitivity analysis
In this work, we performed both local and global sensitivity analyses for the kinetic parameters and initial conditions in the mathematical model. Details of the local sensitivity analysis are presented in Sections 2.5 and 3.2 of the manuscript and the results are depicted in Figure 5 and further Table H: Relative sensitivities determined by local sensitivity analysis of the kinetic parameters  and initial conditions, around the nominal parameter set defined by Tables D and G. The top 10 most sensitive parameters are tabulated for each treatment condition and ordered from most to least sensitive. The relative sensitivity determined by perturbation of +1% is presented in the first column, followed by the relative sensitivity corresponding to a −1% perturbation to the parameter value in the second column, for each treatment condition. Parameters that were also identified as significant by global sensitivity analysis are in boldface (c.f. Figure A).  Table H. Following our previous work [94,95], here we used the Latin hypercube sampling (LHS) method [96] to efficiently sample the parameter space (i.e. to generate the parameter sets and initial conditions) for the global sensitivity analysis. The sampling ranges are specified in Table D for the kinetic parameters related to drug uptake and the regulation of the transcription factors c-Myc and Chop and in Table G for the kinetic parameters related to the intrinsic apoptosis pathway. Using the LHS method, we generated 100,000 sets of parameters/initial conditions and simulated the treatment conditions described in Section 2.2 of the manuscript (untreated, venetoclax monotherapy, tedizolid monotherapy, venetoclax/tedizolid combination treatment) for each set. We used multi-parametric sensitivity analysis (MPSA) [97][98][99], which evaluates the parameter (and initial condition) sensitivities based on Kolmogorov-Smirnov statistics, returning sensitivity values between 0 and 1. A larger parameter sensitivity indicates that the corresponding parameter variation has a large impact on the model output [100]. The results are presented in Figure A for each treatment condition, with the model output taken to be the cell viability at the end of treatment. Inspection of Figure A reveals that the results of the global sensitivity analysis support the results of the local sensitivity analysis results. In particular, kinetic parameters related to the expression of Chop have a significant influence on the response to treatment.